      eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_setup_alpha',num2str(alphaT),...
         '_rho',rhoT,'_approx.mat')]); 
     load MAT/WIOD.mat;
     load('MAT/Step_2_firm_data_alpha4.mat','FI_J','firmsecD_sorted','start_sorted');

% Export exposure
X_mnjf = firmsecD_sorted*(X_mnj(start(france):start(france+1)-1,:));
X_flows = pi_mnjf.*X_mnjf;
turnover_f = sum(X_flows,2);
exposure_X_f = X_flows./repmat(turnover_f,1,N);
exposure_X_f(isnan(exposure_X_f))=0;
% Homogeneous Export exposure
N_jf = sum(firmsecD_sorted,1);
N_jf_repmat = repmat(N_jf,FI_J,1);
one_over_N_repmat = firmsecD_sorted./N_jf_repmat;
one_over_N = sum(one_over_N_repmat,2);
pi_mnjf_hom = repmat(one_over_N,1,N);
X_flows_hom = pi_mnjf_hom.*X_mnjf;
turnover_f_hom = sum(X_flows_hom,2);
exposure_X_f_hom = X_flows_hom ./repmat(turnover_f_hom,1,N);
relative_exposure_X_f = exposure_X_f ./exposure_X_f_hom;
relative_exposure_X_f(exposure_X_f_hom==0)=1;
% corresponds to NT good sectors

% Import exposure
for m=1:N
    exposure_M_f(:,m) = sum(pi_M_f(:,start(m):start(m+1)-1),2);
end
exposure_M_f(exposure_M_f<0)=0;
% Homogeneous Import exposure
one_minus_pi_l_f = 1-pi_l_f;
one_minus_labor_demand_f = one_minus_pi_l_f.*turnover_f;
gamma_pi_X_laborshare_nif = repmat(one_minus_labor_demand_f,[1 J*N]).*pi_M_f;
numerator = firmsecD_sorted'*gamma_pi_X_laborshare_nif;
X_fra = X_mnj(start(france):start(france+1)-1,:);
den = sum(X_fra,2);
denominator = repmat(den,[1 J*N]);
one_minus_alp_gam_hom = numerator./denominator;
clear numerator denominator;
pi_X_laborshare_nif = turnover_f.*pi_l_f;
numerator = firmsecD_sorted'*pi_X_laborshare_nif;
pi_l_hom = numerator./den;
pi_M_hom= one_minus_alp_gam_hom./(1-repmat(pi_l_hom,[1 N*J]));
pi_M_f_hom = firmsecD_sorted*pi_M_hom;
for m=1:N
    exposure_M_f_hom(:,m) = sum(pi_M_f_hom(:,start(m):start(m+1)-1),2);
end
exposure_M_f_hom(exposure_M_f_hom<0)=0;
relative_exposure_M_f = exposure_M_f ./exposure_M_f_hom;
relative_exposure_M_f(exposure_M_f_hom==0)=0;
% Relative size
rho=3*ones(J,1);
rho_f=firmsecD_sorted*rho;
rho_fn = repmat(rho_f,[1 N]);
pi_l_f_repmat = repmat(pi_l_f,[1 N]);
VA_fnj = ((1+pi_l_f_repmat.*(rho_fn-1))./rho_fn).*X_flows;
VA_f = sum(VA_fnj,2);
VA_tot = sum(VA_f);
omega_f_tot = VA_f./VA_tot;
relative_omega_f_tot = omega_f_tot ./(1/FI_J);
repmat_relative_omega_f_tot = repmat(relative_omega_f_tot,1,N);

for m=1:N
    X1 = repmat_relative_omega_f_tot(:,m);
    X2 = relative_exposure_X_f(:,m);
    covX1X2 = cov(X1,X2);
    cov_size_X_relative_exposure(m,1) = covX1X2(1,2);

    X2 = relative_exposure_M_f(:,m);
    covX1X2 = cov(X1,X2);
    cov_size_M_relative_exposure(m,1) = covX1X2(1,2);
end

dlnYn_def = [];
dlnYn_doubledef = [];
dlnYn_homog_def = [];
dlnYn_homog_doubledef = [];
gran_resn_def =[];
gran_resn_doubledef =[];
w_f =[];
dlnY_f = [];

for cty = {'AUS' 'AUT' 'BEL' 'BGR' 'BRA' 'CAN' 'CHN' 'CYP' 'CZE' ...
           'DEU' 'DNK' 'ESP'  'EST' 'FIN' 'GBR' ...
           'GRC' 'HUN' 'IDN' 'IND' 'IRL' 'ITA' 'JPN' 'KOR' ...
           'LTU' 'LVA' 'MEX' 'MLT' 'NLD' 'POL' 'PRT' 'ROU' 'RUS' ...
           'SVK' 'SVN' 'SWE' 'TUR' 'TWN'  'USA' };
    idx = find(ismember(ISO, cty));
    shockT = strcat(cty,'prod_p10');
    shockT=shockT{1};
    
    fn1 = strcat('MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_5_output_alpha',num2str(alphaT),...
          '_rho',rhoT,'_',shockT,'_approx.mat');
    load(fn1,'VA_fj_0','VA_hat_fj','RGDP_def_hat_n','GDP_deflator_n','RGDP_doubledef_hat_n','RGDP_doubledef_hat_fj');
          
    dlnYn_def = [dlnYn_def; (RGDP_def_hat_n(15)-1)];
    dlnYn_doubledef = [dlnYn_doubledef; (RGDP_doubledef_hat_n(15)-1)];
   
    aa = VA_fj_0/sum(VA_fj_0);
    w_f = [w_f aa];
    dlnY_f = [dlnY_f (RGDP_doubledef_hat_fj-1)*100];
    
    clear VA_fj_0 VA_hat_fj RGDP_doubledef_hat_fj RGDP_def_hat_n GDP_deflator_n RGDP_doubledef_hat_n
    
    % Homogeneous
     fn2 = strcat('MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_5_output_alpha',num2str(alphaT),...
           '_rho',rhoT,'_',shockT,'_Hom_approx.mat');
     load(fn2,'RGDP_def_hat_n','RGDP_doubledef_hat_n');   
        
     dlnYn_homog_def = [dlnYn_homog_def; (RGDP_def_hat_n(france)-1)];
     dlnYn_homog_doubledef = [dlnYn_homog_doubledef; (RGDP_doubledef_hat_n(france)-1)];
          
     clear RGDP_doubledef_hat_n RGDP_def_hat_n
     
     fn3 = strcat('MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_6_output_alpha',num2str(alphaT),...
          '_rho',rhoT,'_',shockT,'_approx.mat');
     load(fn3,'cov_e_f_def','cov_e_f_doubledef');
          
     gran_resn_def = [gran_resn_def; cov_e_f_def*.1];
     gran_resn_doubledef = [gran_resn_doubledef; cov_e_f_doubledef*.1];
end

gran_dlnYn_doubledef = gran_resn_doubledef./dlnYn_doubledef;
cov_size_M_relative_exposure = cov_size_M_relative_exposure([1:france-1 france+1:end-1]);
cov_size_X_relative_exposure = cov_size_X_relative_exposure([1:france-1 france+1:end-1]);
cov_size_M_relative_exposure_woGRC = cov_size_M_relative_exposure([1:greece-1 greece+1:end],1); 
cov_size_X_relative_exposure_woGRC = cov_size_X_relative_exposure([1:greece-1 greece+1:end],1); 
gran_dlnYn_doubledef_woGRC = gran_dlnYn_doubledef([1:greece-1 greece+1:end],1); 

labels = {'AUS', 'AUT', 'BEL', 'BGR', 'BRA', 'CAN', 'CHN' ,'CYP', ...
           'CZE', 'DEU', 'DNK', 'ESP' , 'EST', 'FIN',  'GBR', 'GRC', ...
           'HUN', 'IDN', 'IND', 'IRL', 'ITA', 'JPN', 'KOR', ...
           'LTU' ,'LVA', 'MEX','MLT', 'NLD', 'POL', 'PRT', 'ROU','RUS', ...
           'SVK', 'SVN', 'SWE', 'TUR', 'TWN', 'USA'};
figure(8)
plot(cov_size_M_relative_exposure,gran_dlnYn_doubledef,'bo','LineWidth',2)
h = lsline;
set(h(1),'color','r','LineWidth',2)
text(cov_size_M_relative_exposure,gran_dlnYn_doubledef,labels,'VerticalAlignment','bottom','HorizontalAlignment','right')
xlabel('$cov\left(\frac{\omega_{f,n}}{\bar{\omega}},\frac{\pi^M_{f,mn}}{\pi^M_{j,mn}}\right)$','Interpreter','latex')
ylabel('$\Gamma^m/d\ln Y^m$','Interpreter','latex')
fig8 = strcat('Figures/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Fig6a.eps');
saveas(gcf,fig8,'epsc2')

figure(9)
plot(cov_size_X_relative_exposure,gran_dlnYn_doubledef,'bo','LineWidth',2)
h = lsline;
set(h(1),'color','r','LineWidth',2)
text(cov_size_X_relative_exposure,gran_dlnYn_doubledef,labels,'VerticalAlignment','bottom','HorizontalAlignment','right')
xlabel('$cov\left(\frac{\omega_{f,n}}{\bar{\omega}},\frac{s^X_{f,nm}}{s^X_{j,nm}}\right)$','Interpreter','latex')
ylabel('$\Gamma^m/d\ln Y^m$','Interpreter','latex')
fig9 = strcat('Figures/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Fig6b.eps');
saveas(gcf,fig9,'epsc2')
